#Reputation Replication Code
#Ryan Brutger & Josh Kertzer
#Last modified: July 10, 2017

## This R file contains the code necessary to replicate the analysis in the main text. Its companion R file ("Reputation Replication 2.R") replicates the analysis in the supplementary appendix.
# All of the following analyses were carried out using R version 3.4.0 on an 2.9 GHz Intel Core i5 iMac running OS X Sierra version 10.12.5. 

################# 1. Load libraries, load and clean data ##############

#First, set your working directory using the setwd() command:
setwd()

data1 <- get(load("repdata1.RData")) # Data for main analysis
data <- get(load("repdata2.RData")) # Data for STM analysis
data <- subset(repdata2[is.na(repdata2$Casualties)==TRUE | repdata2$Casualties==0, ]) # Subset to hold casualties constant (analysis with alternative casualty levels are discussed in the appendix)

#Load required R packages (download first if necessary)
#Note: the syntax in the tm and stm packages changes periodically; the code below was written to work with tm version 0.7-1 and stm 1.2.2
library(stm)
library(tm) 
library(ggplot2)
library(systemfit) 
library(lmtest)
library(mediation) 

#### Functions for cleaning data
#Rescaling function for more easily interpretable results
rescale <- function(x){
  return((x-min(x,na.rm=TRUE))/(max(x-min(x,na.rm=TRUE),na.rm=TRUE)))
}
#Recoding function
recode <- function(variable, reverse=FALSE, maxVal, binarize=FALSE, x=NA, x2){
	if (is.factor(variable)){variable <- as.numeric(variable)}
	if (missing(maxVal)){maxVal <- max(variable, na.rm=TRUE)}
	variable[variable > maxVal] <- NA
	if (reverse){
		temp <- variable
		for (j in 1:maxVal){
			temp[which(variable==j)] <- maxVal-j+1
		}
		return(temp)
	}
	if (binarize){
		if (is.na(x)){stop("Specify the value to dichotomize on")}
		temp <- variable
		if (missing(x2)){
			i <- which(variable == x)
			k <- which(variable < x | variable > x)
		}
		else{
			i <- which(variable >= x & variable <= x2)
			k <- which(variable < x | variable > x2)
		}
		temp[i] <- 1
		temp[k] <- 0
		return(temp)
	}
	return(variable)
}

#Military assertiveness
#ma1r <- recode(data1$ma1, reverse=TRUE)
#ma3r <- recode(data1$ma3, reverse=TRUE)
#data1$milAssert <- (ma1r + data1$ma2 + ma3r)/3 # Variable already saved in dataset, coding scale shown here as demonstration

#National chauvinism
#nc1r <- recode(data1$nc1, reverse=TRUE)
# data1$natChauv <- (nc1r + data$nc2)/2  # Variable already saved in dataset, coding scale shown here as demonstration

#Rescale and code main variables of interest
data1$milAssert1 <- rescale(data1$milAssert)
data1$natChauv1 <- rescale(data1$natChauv)
data1$intTrust1 <- rescale(data1$intTrust)
data1$Republican <- ifelse(data1$polParty==1,1,ifelse(data1$polParty>1,0,NA))
data1$Democrat <- ifelse(data1$polParty==2,1,ifelse(data1$polParty==1,0,ifelse(data1$polParty>2,0,NA)))

##Code new education variable to a four-level ordinal variable (high school or less, some college, college/university, postgraduate)
data1$education2 <- ifelse(data1$education==1| data1$education==4 | data1$education==10, 1, 
                           ifelse(data1$education==3 | data1$education==5 | data1$education==7, 2,
                                  ifelse(data1$education==2, 3,
                                         ifelse(data1$education==6 | data1$education==8, 4, NA))) )

data1$Male <- ifelse(data1$gender==1,1,ifelse(data1$gender==2,0,NA))

################# 2. Main experimental results (ATEs and regression models) ##############

set.seed(43215)

#Begin by looking for bootstrapped treatment effects
#Function for calculating bootstrapped quantities of interest
bootTreat2 <- function(dat, dv="totApprov", B=1500){
  bootResults <- matrix(NA, nrow=B, ncol=5)
  for (i in seq_len(B)){
    resample <- sample(1:nrow(dat),nrow(dat),replace=T)
    temp <- dat[resample,]
    dep.var <- eval(parse(text=paste("temp",dv,sep="$")))
    A <- mean(dep.var[which(temp$notEngage3==1)], na.rm=TRUE)
    B <- mean(dep.var[which(temp$engage3==1)], na.rm=TRUE)
    C <- mean(dep.var[which(temp$stayOut3==1)], na.rm=TRUE)
    bootResults[i,1] <- (A-C)
    bootResults[i,2] <- (B-C)
    bootResults[i,3] <- (A-B)
    bootResults[i,4] <- (B-C)/(A-C)
    bootResults[i,5] <- (A-B)/(A-C)
    drop(list())
  }
  return(list(dv=dv, boot=bootResults, Audience.Cost=mean(bootResults[,1], na.rm=TRUE), Force.Cost=mean(bootResults[,2], na.rm=TRUE), Inconsistency.Cost=mean(bootResults[,3], na.rm=TRUE), Force.Fraction=mean(bootResults[,4], na.rm=TRUE), Inconsistency.Fraction=mean(bootResults[,5], na.rm=TRUE)))
}

## Look at ATE for reputation costs
## Country-level reputation costs
arep.full <- bootTreat2(data1, dv="AmerRep") 
arep.full$Inconsistency.Cost #Inconsistency cost: 0.711 (18% shift)
arep.full$Force.Cost #Belligerence cost: -0.029
#To calculate p-values:
1 - length(which(arep.full$boot[,3] > 0))/nrow(arep.full$boot) #p < 0.000 for inconsistency reputation cost
1 - length(which(arep.full$boot[,2] < 0))/nrow(arep.full$boot) #p < 0.394 for belligerence reputation cost 

## President-level reputation costs
prep.full <- bootTreat2(data1, dv="PresRep")
prep.full$Inconsistency.Cost #Inconsistency cost: 0.755 (19% shift)
prep.full$Force.Cost #Belligerence cost: -0.099
#To calculate p-values:
1 - length(which(prep.full$boot[,3] > 0))/nrow(prep.full$boot) #p < 0.00 for inconsistency reputation cost
1 - length(which(prep.full$boot[,2] < 0))/nrow(prep.full$boot) #p < 0.188 for belligerence reputation cost

###### Table 1: Seemingly Unrelated Regressions (SURs) ########

## Check correlation between country and president's reputation (see appendix for further discussion and analysis)
round(cor(data1$AmerRep, data1$PresRep, use="complete.obs"),digits=2) ## r = 0.86

dat.sys <- data1[which(!is.na(data1$AmerRep) & !is.na(data1$PresRep)),] #First create common dataframe

amod.1 <- AmerRep ~ stayOut3 + notEngage3 + milAssert1 + intTrust1 + natChauv1 + Democrat + Republican + education2 + income + Male
pmod.1 <- PresRep ~ stayOut3 + notEngage3 + milAssert1 + intTrust1 + natChauv1 + Democrat + Republican + education2 + income + Male
sys.1 <- systemfit(method="SUR", list(amod.1, pmod.1), data=dat.sys)
summary(sys.1) #N: 586, Adj. R2: 0.119; #N: 586, Adj. R2: 0.127

#Produce output table (slightly ugly)
l <- length(coef(sys.1)) #22
coefA <- round(coef(sys.1)[1:(l/2)], digits=3)
coefP <- round(coef(sys.1)[((l/2)+1):l], digits=3)
seA <- round(sqrt(diag(vcov(sys.1)[1:(l/2),1:(l/2)])), digits=3)
seP <- round(sqrt(diag(vcov(sys.1)[((l/2)+1):l,((l/2)+1):l])), digits=3)
#Add parentheses
seA.2 <- paste("(",seA,")",sep="")
seP.2 <- paste("(",seP,")",sep="")
#Since belligerence cost is actually engage - stay out, rather than stay out - engage, flip the sign of the second coefficient in each model
coefA[2] <- coefA[2]*-1
coefP[2] <- coefP[2]*-1

outMat <- cbind(as.vector(rbind(coefA, seA.2)), as.vector(rbind(coefP, seP.2)))
rownames(outMat) <- rep("",22)
rownames(outMat)[seq(1,22,by=2)] <- c("Intercept", "Belligerence", "Inconsistency", "Militant assertiveness", "International trust", "National chauvinism", "Democrat", "Republican", "Education", "Income", "Male")
colnames(outMat) <- c("Country", "President")
outMat #Output table

#Do wald test to compare fits (footnote 72)
#We need dataframes of the same size
wald.dat <- data1[which(!is.na(data1$AmerRep) & !is.na(data1$PresRep) & !is.na(data1$milAssert1) & !is.na(data1$intTrust1) & !is.na(data1$natChauv1) & !is.na(data1$Democrat) & !is.na(data1$Republican) & !is.na(data1$education2) & !is.na(data1$income) & !is.na(data1$Male)),]
amod.1.lm <- lm(AmerRep ~ stayOut3 + notEngage3 + milAssert1 + intTrust1 + natChauv1 + Democrat + Republican + education2 + income + Male, data=wald.dat)
amod.2.lm <- lm(AmerRep ~ stayOut3 + notEngage3, data=wald.dat)
waldtest(amod.1.lm, amod.2.lm) #F = 3.792, p < 0.000

pmod.1.lm <- lm(PresRep ~ stayOut3 + notEngage3 + milAssert1 + intTrust1 + natChauv1 + Democrat + Republican + education2 + income + Male, data=wald.dat)
pmod.2.lm <- lm(PresRep ~ stayOut3 + notEngage3, data=wald.dat)
waldtest(pmod.1.lm, pmod.2.lm) #F = 4.213, p < 0.000

######### Figure 1: Hawks and doves assess reputation costs in significantly different ways #########

set.seed(43215)

## Split sample into hawks and doves based on military assertiveness quartiles
quantile(data1$milAssert1, c(0.25, 0.75), na.rm=TRUE) #0.42, 0.67
data1$milAssert2 <- ifelse(data1$milAssert1 <= 0.42, 0,ifelse(data1$milAssert1 >= 0.66,1,NA))

## America's Reputation
## Doves
data1.doves <- subset(data1[data1$milAssert2==0 & is.na(data1$milAssert2)==FALSE, ])
full.doves <- bootTreat2(data1.doves, dv="AmerRep") 

## Hawks
data1.hawks <- subset(data1[data1$milAssert2==1 & is.na(data1$milAssert2)==FALSE, ])
full.hawks <- bootTreat2(data1.hawks, dv="AmerRep") 

## President's Reputation
## Doves

full.doves.prep <- bootTreat2(data1.doves, dv="PresRep") 

## Hawks
full.hawks.prep <- bootTreat2(data1.hawks, dv="PresRep") 

## Density plots for Hawks versus Doves, currently reporting fraction of total reputation cost

B <- 1500 #Number of bootstraps

dens <- data.frame(y=c(arep.full$boot[,2], arep.full$boot[,3], full.doves$boot[,2], full.doves$boot[,3],  full.hawks$boot[,2], full.hawks$boot[,3], prep.full$boot[,2], prep.full$boot[,3], full.doves.prep$boot[,2], full.doves.prep$boot[,3], full.hawks.prep$boot[,2], full.hawks.prep$boot[,3]), dv=rep(c("America's Reputation", "President's Reputation"), each=6*B), samp=rep(rep(rep(c("Full Sample", "Doves", "Hawks"), each=B),2),each=2), Cost=rep(rep(c("Belligerence", "Inconsistency"), each=B),6))  

#Reorder levels
dens$samp <- factor(dens$samp, levels=c("Full Sample", "Doves", "Hawks"))

#Create density plot
dev.new(height=4.5, width=11)
ggplot(dens, aes(x=y)) + geom_density(aes(fill=Cost), alpha=0.35) + facet_grid(dv ~ samp, scales="fixed") + labs(x="Reputation cost treatment effect", y="Density") + theme_bw() + theme(axis.text.y=element_blank(), axis.ticks.y=element_blank()) + geom_vline(xintercept=0, linetype=3, alpha=0.3) 
#Because ggplot is sometimes finicky about color settings, you can convert to black and white in any image-editing program (e.g. Adobe Illustrator). To increase legibility for the printed version, we convert the colors to RGB (76,76,76) and (232,232,232), and change the opacity to 90% and 50%, respectively.

############ 3. Text analysis using structural topic models (STMs) ###########

##### Estimating STMs necessary to produce Figures 2 and 3 ######

## Clean and preprocess the text before estimating the STM; if this produces errors, be sure to check what version of the tm and stm package you have installed
new.data <- data[c("Territory",   "Strategy", "totApprov", "AmerRep", "PresRep", "milAssert", "amerRepFree")] # America's reputation
new.data2 <- data[c("Territory",   "Strategy", "totApprov", "AmerRep", "PresRep", "milAssert", "presRepFree")] # President's reputation
meta.data <- data[c("Territory",   "Strategy", "totApprov", "AmerRep", "PresRep", "milAssert")]
docs1 <- data[c("amerRepFree")]
docs2 <- data[c("presRepFree")]
Corpus <- Corpus(VectorSource(docs1$amerRepFree))
Corpus2 <- Corpus(VectorSource(docs2$presRepFree))

## Clean corpus
Corpus.clean <- tm_map(Corpus, stripWhitespace)
Corpus.clean <- tm_map(Corpus.clean, removePunctuation)
Corpus.clean <- tm_map(Corpus.clean, tolower)

Corpus.clean2 <- tm_map(Corpus2, stripWhitespace)
Corpus.clean2 <- tm_map(Corpus.clean2, removePunctuation)
Corpus.clean2 <- tm_map(Corpus.clean2, tolower)

dataframe<-data.frame(text=unlist(sapply(Corpus.clean, `[`)), stringsAsFactors=F)
dataframe2<-data.frame(text=unlist(sapply(Corpus.clean2, `[`)), stringsAsFactors=F)

clean.data <- cbind(meta.data, dataframe)
clean.data2 <- cbind(meta.data, dataframe2)
colnames(clean.data) <- c("Territory",   "Strategy", "totApprov", "AmerRep", "presRep", "milAssert","documents")
colnames(clean.data2) <- c("Territory",   "Strategy", "totApprov", "AmerRep", "presRep", "milAssert","documents")

## combine data so all reputation responses are included 
clean.data3 <- rbind(clean.data, clean.data2)

## Remove offending characters that cause errors
for(i in 1:nrow(clean.data3)){
  clean.data3$documents[i] <- gsub("[^[:alnum:]///' ]", "", clean.data3$documents[i])
}

## Now create dataframes with just doves or hawks for each treatment group
quantile(clean.data3$milAssert, c(0.25, 0.75), na.rm=TRUE) ## Quartile cutpoints 2.33 and 3.33
clean.data3$milAssert1 <- ifelse(clean.data3$milAssert <= 2.34, 0,ifelse(clean.data3$milAssert >= 3.33,1,NA))

clean.data4 <- clean.data3[c("Territory",   "Strategy", "totApprov", "AmerRep", "presRep", "milAssert1", "documents")] 
clean.data4 <- subset(clean.data4[is.na(clean.data4$milAssert1)==FALSE, ])

## Use "DH" to identify data for comparing doves and hawks
DH.stayout <- subset(clean.data4[clean.data4$Strategy==1, ])
DH.engage <- subset(clean.data4[clean.data4$Strategy==2, ])
DH.notengage <- subset(clean.data4[clean.data4$Strategy==3, ])

## Doves versus Hawks in StayOut
processed <- textProcessor(DH.stayout$documents, metadata=DH.stayout)
out <- prepDocuments(processed$documents, processed$vocab, processed$meta)

## Now run manyTopics to aid in model selection
# output is used to select the number of topics, but the output is not reported in the paper or appendix; you can also skip right to line 285
DH.stayout.multi <- manyTopics(processed$documents, processed$vocab, K = 6:16, prevalence = ~ milAssert1, max.em.its = 300, runs = 20, data = processed$meta, seed = 51111)

t <- DH.stayout.multi$semcoh
s <- DH.stayout.multi$exclusivity
sc <- c()
ex <- c()
cols <- c()
for(i in 1:11){
  sc <- c(sc, mean(t[[i]]))
  ex <- c(ex, mean(s[[i]]))
}
plot(sc, ex, xlim = c(-150, -205),
     ylim = c(9.0, 10.0), color = cols)

ex # 9.206221 9.279210 9.380158 9.434291 9.503567 9.552080 9.566482 9.631469 9.649919 9.674960 9.662049
# Exclusivity plateaus around 9.6 once 13 topics are used. We choose to analyze 13 topics since it achieves high exclusivity, which is helpful for observing differences in topics, and relatively little exclusivity is gained by adding additional topics. Although we sacrifice semantic coeherence by focusing on exclusivity, we still only find only one topic with a significant difference that lacks a clear substantive interpretation. 

## Doves versus Hawks in Engage (pick up here if you're skipping the manyTopics function)
processed <- textProcessor(DH.engage$documents, metadata=DH.engage)

out <- prepDocuments(processed$documents, processed$vocab, processed$meta)

## Fit the structural topic model
DH.engage.model <- selectModel(processed$documents, processed$vocab, K=13, 
                               prevalence = ~ milAssert1, max.em.its = 300,
                               runs=20, data = processed$meta, seed = 51111)


DH.engage.model1 <- DH.engage.model$runout[[2]] #Select second model run

plot(DH.engage.model1, type="summary", xlim=c(0,1), n=7)

DH.engage.model1.effect <-  estimateEffect(1:13 ~ milAssert1, DH.engage.model1, meta=processed$meta, uncertainty = "Global")

####### Figure 2: Difference in Topic Proportion Between Doves and Hawks: #####

### Left-hand panel: Effect of Doves vs Hawks in Engage 
plot(DH.engage.model1.effect, covariate = "milAssert1", topics=c(6, 11, 10), method = "difference", cov.value1=1, cov.value2=0, xlab="Change in topical prevalence from Dove to Hawk", main="Effect of Doves vs Hawks in Engage", xlim=c(-.1,.1),  labeltype = "custom", custom.labels = c('Accepting Unpopularity', 'Anti-Interventionism', 'Presidential Deference'))
#(topic labels are moved in the paper to accomodate margins)

### To speed up estimation, the code for the right-hand panel of Figure 2 is presented beginning on line 328

######### Figure 3  ###########
## Representative responses presented in the left-hand column ("Engage") of Figure 3. For the right-hand column, see code beginning on line 333 
#The figure was manually assembled outside of R and responses were reconstituted to their original form prior to cleaning for STM
## The following idendifies 10 representative responses for each topic, from which the authors selected the responses displayed in the paper

#Accepting Unpopularity 
engage.thoughts6 <- findThoughts(DH.engage.model1, texts=processed$meta$documents, topics=6, n=10)
engage.thoughts6
#Anti-Interventionism
engage.thoughts11 <- findThoughts(DH.engage.model1, texts=processed$meta$documents, topics=11, n=10)
engage.thoughts11
#Presidential Deference
engage.thoughts10 <- findThoughts(DH.engage.model1, texts=processed$meta$documents, topics=10, n=10)
engage.thoughts10

### Now, results for the Not Engage condition

## Doves versus Hawks in NotEngage
processed <- textProcessor(DH.notengage$documents, metadata=DH.notengage)

out <- prepDocuments(processed$documents, processed$vocab, processed$meta)

## Fit the structural topic model
DH.notengage.model <- selectModel(processed$documents, processed$vocab, K=13, 
                                  prevalence = ~ milAssert1, max.em.its = 300,
                                  runs=20, data = processed$meta, seed = 51111)

DH.notengage.model1 <- DH.notengage.model$runout[[1]]

set.seed(43215)
DH.notengage.model1.effect <-  estimateEffect(1:13 ~ milAssert1, DH.notengage.model1, meta=processed$meta, uncertainty = "Global")


####### Figure 2: Difference in Topic Proportion Between Doves and Hawks: #####

### Right-hand panel: Effect of Doves vs Hawks in Not Engage 
plot(DH.notengage.model1.effect, covariate = "milAssert1", topics=c(4, 5, 11), method = "difference", cov.value1=1, cov.value2=0, xlab="Difference of change from Dove to Hawk", main="Effect of Doves vs Hawks in Not Engage", xlim=c(-.1,.1),  labeltype = "custom", custom.labels = c('Accepting Unpopularity', 'Anti-Interventionism', 'Inconsistency')) #Topic labels are moved in the paper to accomodate margins

######### Figure 3  ###########

## Representative responses presented in the right-hand column ("Not Engage") of Figure 3.  
#The Figure was manually assembled outside of R and responses were reconstituted to their original form (prior to cleaning for STM)
## The following identifies 10 representative responses for each topic, from which the authors selected the responses displayed in the paper

#Accepting Unpopularity 
notengage.thoughts4 <- findThoughts(DH.notengage.model1, texts=processed$meta$documents, topics=4, n=10)
notengage.thoughts4
#Anti-Interventionism
notengage.thoughts5 <- findThoughts(DH.notengage.model1, texts=processed$meta$documents, topics=5, n=10)
notengage.thoughts5
#Presidential Deference
notengage.thoughts11 <- findThoughts(DH.notengage.model1, texts=processed$meta$documents, topics=11, n=10)
notengage.thoughts11

## Calculate estimated topic proportions for Figure 3

apply(make.dt(DH.engage.model1),2,mean)[1+c(6,11,10)] #Extract just results for topics 6,11 and 10
#Topic 6 (accepting unpopularity): 8%
#Topic 11 (anti-interventionism): 8%
#Topic 10 (presidential deference): 7%

apply(make.dt(DH.notengage.model1),2,mean)[1+c(4,5,11)] #Extract just results for topics 4,5 and 11

#Topic 4 (accepting unpopularity): 5%
#Topic 5 (anti-interventionism): 8%
#Topic 11 (inconsistency): 10%

################# 4. Moderated mediation analysis  ###################

########### Figure 4: Reputational Consequences of Belligerence and Inconsistency on Presidential Approval #######

set.seed(43215)

#Isolate belligerence and inconsistency treatments, and subset based on hawks versus doves
dat.bel <- subset(data1, data1$notEngage3==0)
dat.inc <- subset(data1, data1$stayOut3==0)
dat.bel.dove <- subset(dat.bel, dat.bel$milAssert2==0)
dat.bel.hawk <- subset(dat.bel, dat.bel$milAssert2==1)
dat.inc.dove <- subset(dat.inc, dat.inc$milAssert2==0)
dat.inc.hawk <- subset(dat.inc, dat.inc$milAssert2==1)

## America's Reputation - Belligerence - Doves
arep.med <- lm(AmerRep ~ engage3 + education2 + income + Male + ideology, data=dat.bel.dove)
arep.dv <- lm(totApprov ~ AmerRep + engage3 + education2 + income + Male + ideology, data=dat.bel.dove)
arep.med.b.d <- mediate(arep.med, arep.dv, treat="engage3", mediator="AmerRep", sims=1500)
summary(arep.med.b.d) #ACME: -0.4081, p < 0.02

## America's Reputation - Belligerence - Hawks
arep.med <- lm(AmerRep ~ engage3 + education2 + income + Male + ideology, data=dat.bel.hawk)
arep.dv <- lm(totApprov ~ AmerRep + engage3 + education2 + income + Male + ideology, data=dat.bel.hawk)
arep.med.b.h <- mediate(arep.med, arep.dv, treat="engage3", mediator="AmerRep", sims=1500)
summary(arep.med.b.h) #ACME: 0.3964, p < 0.08

## America's Reputation - Inconsistency - Doves
arep.med <- lm(AmerRep ~ notEngage3 + education2 + income + Male + ideology, data=dat.inc.dove)
arep.dv <- lm(totApprov ~ AmerRep + notEngage3 + education2 + income + Male + ideology, data=dat.inc.dove)
arep.med.i.d <- mediate(arep.med, arep.dv, treat="notEngage3", mediator="AmerRep", sims=1500)
summary(arep.med.i.d) #ACME: -0.393, p < 0.11

## America's Reputation - Inconsistency - Hawks
arep.med <- lm(AmerRep ~ notEngage3 + education2 + income + Male + ideology, data=dat.inc.hawk)
arep.dv <- lm(totApprov ~ AmerRep + notEngage3 + education2 + income + Male + ideology, data=dat.inc.hawk)
arep.med.i.h <- mediate(arep.med, arep.dv, treat="notEngage3", mediator="AmerRep", sims=1500)
summary(arep.med.i.h) #ACME: -1.683, p < 0.01

# Belligerence values - America's Reputation
points.b.d <- c(arep.med.b.d$tau.coef , arep.med.b.d$z.avg, arep.med.b.d$d.avg)
b.d.CI.low <- c(arep.med.b.d$tau.ci[1],arep.med.b.d$z.avg.ci[1],arep.med.b.d$d.avg.ci[1])
b.d.CI.high <- c(arep.med.b.d$tau.ci[2],arep.med.b.d$z.avg.ci[2],arep.med.b.d$d.avg.ci[2]) 
points.b.h <- c(arep.med.b.h$tau.coef , arep.med.b.h$z.avg, arep.med.b.h$d.avg)
b.h.CI.low <- c(arep.med.b.h$tau.ci[1], arep.med.b.h$z.avg.ci[1], arep.med.b.h$d.avg.ci[1])
b.h.CI.high <- c(arep.med.b.h$tau.ci[2],arep.med.b.h$z.avg.ci[2], arep.med.b.h$d.avg.ci[2]) 
labels <- c("Total Effect", "ADE", "ACME")

# Inconsistency values - America's Reputation
points.i.d <- c(arep.med.i.d$tau.coef , arep.med.i.d$z.avg, arep.med.i.d$d.avg)
i.d.CI.low <- c(arep.med.i.d$tau.ci[1],arep.med.i.d$z.avg.ci[1],arep.med.i.d$d.avg.ci[1])
i.d.CI.high <- c(arep.med.i.d$tau.ci[2],arep.med.i.d$z.avg.ci[2],arep.med.i.d$d.avg.ci[2]) 
points.i.h <- c(arep.med.i.h$tau.coef , arep.med.i.h$z.avg, arep.med.i.h$d.avg)
i.h.CI.low <- c(arep.med.i.h$tau.ci[1], arep.med.i.h$z.avg.ci[1], arep.med.i.h$d.avg.ci[1])
i.h.CI.high <- c(arep.med.i.h$tau.ci[2],arep.med.i.h$z.avg.ci[2], arep.med.i.h$d.avg.ci[2]) 
labels <- c("Total Effect", "ADE", "ACME")

#Combine quantities of interest
temp <- data.frame(X=c(points.b.d, points.b.h, points.i.d, points.i.h), Low=c(b.d.CI.low, b.h.CI.low, i.d.CI.low, i.h.CI.low), High=c(b.d.CI.high, b.h.CI.high, i.d.CI.high, i.h.CI.high), Cost=(rep(c("Belligerence", "Inconsistency"), each=length(points.b.d)*2)), Group=rep(c("Doves", "Hawks"), each=3), Z=(rep(c("Total", "ADE", "ACME"),4)))

#Now generate plot
dev.new(height=7, width=7)
ggplot(data=temp, aes(x=X, y=Group)) + geom_point() + geom_segment(aes(x=Low, y=Group, xend=High, yend=Group)) + facet_grid(Z~Cost) + theme_bw() + geom_line(aes(x=X,y=Group, group=Cost, alpha=0.8)) + labs(x="Effect", y="") + geom_vline(xintercept=0, linetype=3) + guides(alpha=FALSE)

##### Now, calculate formal moderated mediation p-values to manually annotate onto plot
set.seed(43215)

#Belligerence costs, America's reputation
arep.mod.med <- lm(AmerRep ~ engage3*milAssert2 + education2 + income + Male + ideology, data=dat.bel)
arep.mod.dv <- lm(totApprov ~ engage3*milAssert2 + AmerRep + education2 + income + Male + ideology, data=dat.bel)
mod.med <- mediate(arep.mod.med, arep.mod.dv, sims=1000, boot=FALSE, treat="engage3", mediator="AmerRep", data=dat.bel)
modmed <- test.modmed(mod.med, covariates.1=list(milAssert2=1), covariates.2=list(milAssert2=0), sims=1500)
modmed # ACME(hawks) - ACME(doves) = 1.02, p < 0.001 for ACME; ADE(hawks) - ADE(doves) = 1.00, p < 0.049 for ADE

#Inconsistency costs, America's reputation
arep.mod.med <- lm(AmerRep ~ notEngage3*milAssert2 + education2 + income + Male + ideology, data=dat.inc)
arep.mod.dv <- lm(totApprov ~ notEngage3*milAssert2 + AmerRep + education2 + income + Male + ideology, data=dat.inc)
mod.med <- mediate(arep.mod.med, arep.mod.dv, sims=1000, boot=FALSE, treat="notEngage3", mediator="AmerRep", data=dat.inc)
modmed <- test.modmed(mod.med, covariates.1=list(milAssert2=1), covariates.2=list(milAssert2=0), sims=1500)
modmed # ACME(hawks) - ACME(doves) = -0.710, p < 0.057 for ACME; ADE(hawks) - ADE(doves) = -0.683, p < 0.141 for ADE

#Also superimposed onto the plot: tests of whether the total effects significantly differ from each other (equivalent to the p-value of the interaction):

## America's reputation
#Belligerence
arep.mod.tot <- lm(totApprov ~ engage3*milAssert2 + education2 + income + Male + ideology, data=dat.bel)
summary(arep.mod.tot) #p < 0.001

#Inconsistency
arep.mod.tot <- lm(totApprov ~ notEngage3*milAssert2 + education2 + income + Male + ideology, data=dat.inc)
summary(arep.mod.tot) #p < 0.021

##
## Footnotes 85-86: Now calculate difference in ACMEs using the continious measure of militant assertiveness
##

set.seed(43215)

#Belligerence costs, America's reputation - Continous measure of military assertiveness
arep.mod.med <- lm(AmerRep ~ engage3*milAssert1 + education2 + income + Male + ideology, data=dat.bel)
arep.mod.dv <- lm(totApprov ~ engage3*milAssert1 + AmerRep + education2 + income + Male + ideology, data=dat.bel)
mod.med <- mediate(arep.mod.med, arep.mod.dv, sims=1000, boot=FALSE, treat="engage3", mediator="AmerRep", data=dat.bel)
modmed <- test.modmed(mod.med, covariates.1=list(milAssert1=1), covariates.2=list(milAssert1=0), sims=1000)
modmed # 1.582, p < 0.001 for ACME; 1.990, p < 0.008 for ADE

#Inconsistency costs, America's reputation - Continous measure of military assertiveness
arep.mod.med <- lm(AmerRep ~ notEngage3*milAssert1 + education2 + income + Male + ideology, data=dat.inc)
arep.mod.dv <- lm(totApprov ~ notEngage3*milAssert1 + AmerRep + education2 + income + Male + ideology, data=dat.inc)
mod.med <- mediate(arep.mod.med, arep.mod.dv, sims=1000, boot=FALSE, treat="notEngage3", mediator="AmerRep", data=dat.inc)
modmed <- test.modmed(mod.med, covariates.1=list(milAssert1=1), covariates.2=list(milAssert1=0), sims=1000)
modmed # -1.576, p < 0.002 for ACME; -1.596, p < 0.014 for ADE
